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On the basis of pertubative QCD and the relativistic quark model we calculate 
relativistic corrections to the process of pair J/ip production in proton-proton colli- 
sions at LHC energy yS = 7 TeV. Relativistic terms in the production amplitude 
connected with the relative motion of heavy quarks and the transformation law of 
the bound state wave functions to the reference frame of moving J /ip mesons are 
taken into account. For the gluon and quark propagators entering the amplitude 
we use a truncated expansion in relative quark momenta up to the second order. 
Relativistic corrections to the quark bound state wave functions are considered by 
means of the Breit-like potential. It turns out that the examined effects decrease 
initial nonrelativistic cross section more than two times. The final result lies below 
the experimental value measured by LHCb. 
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I. INTRODUCTION 

The production of heavy quarkonium states in different reactions is the subject of con- 
siderable interest during last years. The mechanism of heavy quarkonium production repre- 
sents the long-standing problem of quantum chromodynamics 0-0|- Most current theoret- 
ical investigations are performed on the basis of nonrelativistic quantum chromodynamics 
(NRQCD) |5| and the quark models. According to these approaches the production of heavy 
quarkonium is divided into two stages. On the first stage one or several quark- ant iquark 
pairs are produced at small distances of order 1/rriQ. This short-range part in the produc- 
tion is associated with the basic interaction of free quarks and gluons and can be evaluated 
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by perturbation theory. The subsequent nonperturbative transition from the intermediate 
states of quarks QQ... and antiquarks QQ... to physical quarkonium states, on the second 
stage, involves long-distance scales of order of quarkonium size l/(m,Qv). The formation of 
the quark bound states is parameterized by nonperturbative matrix elements in NRQCD 
and calculated by means of the bound state wave functions in the quark models. The finding 
of the correspondence between parameters of the quark models and NRQCD opens the way 
for better understanding of the quark-gluon dynamics. Both approaches complement each 
other and can reveal new aspects of the color dynamics of quarks and gluons. 

One of the directions in this investigation is related with the pair production of double 
heavy quarkonium. The initial impulse to intensive investigations was given in this field 
several years ago by the measurement of the pair charmonium production cross sections 
in electron-positron annihilation. The experimental data obtained at the Belle and BaBar 
detectors disagreed with the calculations based on NRQCD. The theoretical results were 
improved and adjusted in correspondence with the experiment after the account of one-loo 
perturbative corrections and relativistic corrections to the nonrelativistic cross section 
|8[. One of the learned lessons from this problem consists in the understanding that only 
sequential relativistic approach to the heavy quarkonium production processes can lead 
to reliable theoretical results. It is necessary to point out that subsequent observation 
of numerous charmonium-like states by the Belle and BaBar collaborations with unclear 
nature demands further continuation of the investigations in this direction j9[ . Recently, the 
first experimental result of the LHCb collaboration on the pair charmonium production in 
proton-proton interaction was published jl 



«c b = 5.1±1.0±l.lnb, (1) 

where the first uncertainty is statistical and the second systematic. It agrees with the 
theoretical estimation of the cross section o = 17 4- 22 nb obtained in the leading order of 
QCD where the process of the gluon fusion is the dominant one (ll-15||. These calculations 



give the total value of the cross section for the pair charmonium production a = 3 -j- 5 nb 
in the kinematical region of the LHCb experiment (the region of rapidities 2 < y < 4.5). 
The theoretical uncertainty remains sufficiently large. To the present the calculations of the 
pair charmonium production in pp interaction were carried out in the leading order of QCD 
without inclusion of relativistic corrections. In addition to these permanent theoretical errors 
known from the experience of e + e~ annihilation we have in this task new specific uncertainty 
caused by parton distribution functions at small x values because the gluon contribution 
from the region of small x is dominant. There appears also another uncertainty related with 
the double parton interaction [16]. In this work we study one aspect of the improvement 
of the previous calculations connected with the account of relativistic corrections. Using 



the methods of the relativistic quark model [g, |17H19| we perform new calculation of the 
cross section a{pp —> J/ipJ/ip) accounting for relativistic corrections to the production 
amplitude and the bound state wave functions of heavy c-quarks. So, the aim of this study 
consists in the relativistic description of the pair charmonium production at hadron colliders. 
It is important to note that the interest to the inclusive reactions p + p — > 2J/ip + X, 
pN — > 2J/ip + X is not limited only by the investigation of the production mechanism. The 
study of the quarkonium production in the nuclear matter leads to new data about QCD at 
high density and temperature 0, 2l| . 



There exist different mechanisms for the pair charmonium production in pp-collisions. 
At the collider energies, double quarkonium production occurs through the gluon-gluon 
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interaction channel. In the color singlet model (CSM) a pair of quark-antiquark (cc) is 
created at short distances in color singlet state. Then it evolves nonperturbatively into an 
observed meson J /if}. At small transverse momenta and small invariant masses of the J /if) 
pair the color singlet mechanism gives the main contribution to the cross section. Another 
possibility is to create a pair (cc) with different color and angular momentum from that 
of the final meson. Then the color octet pair evolves to the color singlet charmonium 
emitting soft gluons. This color octet mechanism (COM) plays significant role in the region 
of high transverse momentum. Among large number of the Feynman diagrams describing 
the production of a pair J / if) it has been separated a class in which the J / if) pair production 
is related with the double gluon fragmentation. In this study we analyze the total set of the 
production amplitudes in the color singlet model making primary emphasis upon relativistic 
effects. 



II. GENERAL FORMALISM 

The differential cross section da for the inclusive double charmonium production in 
proton-proton interaction can be presented in the form of the convolution of partonic cross 
section da(gg — > J /if) J /if)) with the parton distribution functions (PDF) in the initial pro- 
tons [ll,[2|: 

da[p + p -»■ J/ip + J/ip + X] = J dxtdx 2 f g/p (xi, ii)f g/p (x 2 , n) da[gg J/if)J/ip], (2) 

where f g / p (x, //) is partonic distribution function for the gluon in the proton, x± )2 are longi- 
tudinal momentum fractions of gluons, /i is the factorization scale. Neglecting the proton 
mass and taking the cm. reference frame of initial protons with the beam along the z- 
axis we can present the gluon on mass-shell momenta fci j2 = £1,2-^(1, 0, 0, ±1). y/S is the 
center-of-mass energy in proton-proton collision. 

In the quasipotentional approach the invariant transition amplitude for the gluonic sub- 
process g + g —7- J /if) + J /if) can be presented as a convolution of a perturbative production 
amplitude of two c-quark and c-antiquark pairs T(pi,p 2 ; qi, 92) and the quasipotential wave 
functions of the final mesons ^ j/^ |8|: 

M[gg^J/ipj/if)](k u k 2 ,P,Q) = J J P)${q, Q) (g) T{pi,p 2 ; gi, q 2 ), (3) 

where p\ and p 2 are four-momenta of c-quark and c-antiquark in the pair forming the first 
J/ip particle, and q 2 and q\ are appropriate momenta for quark and antiquark in the second 
meson J /if). They are defined in terms of total momenta P(Q) and relative momenta p(q) 
as follows: 

Pi,2 = \p±P, (pP)=0; q 1>2 = ±Q±q, {qQ) = 0, (4) 

In Eq. ((3]) we integrate over the relative three-momenta of quarks and antiquarks in the 
final state. The systematic account of all terms depending on the relative quark momenta 
p and q in ([2]) is important for increasing the accuracy of the calculation, p = Lp(0, p) and 
q = Lq(0, q) are the relative four-momenta obtained by the Lorentz transformation of four- 
vectors (0, p) and (0, q) to the reference frames moving with the four-momenta P and Q. 
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FIG. 1: The typical LO diagrams contributing to the partonic process g + g — > J/ip + J /if)- The 
others can be obtained by reversing the quark lines or interchanging the initial gluons. 



The parton- level differential cross section for g + g — >■ J/if)+J/if> is expressed further through 
the Mandelstam variables s, t and u: 

s= (h + k 2 ) 2 = (P + Q) 2 = Xl x 2 S, (5) 

t = (P-h) 2 = (Q-k 2 ) 2 = M 2 -x l VS(P -\P\cos(j)) = M 2 -x 1 x 2 S+x 2 VS(P +\P\cos(f)), 
u = (P-k 2 ) 2 = (Q-h) 2 = M 2 -x 2 y/S(P + \P\cos<f)) = M 2 -x 1 x 2 S+x 1 \ r S(P -\~P\cos(j)), 

where M is the charmonium mass, (f) is the angle between P and the z-axis. The transverse 
momentum Pt of J/ip and its energy Pq can be written as 

p2 I-DI2 • 2, . (M 2 ~t) 2 X ± X 2 \^S X X -X 2 , , . v 

P£ — P sin <p — — t — , Pq — 1 |P|cos0. (6) 

X\X 2 S Xi + X 2 X\ + x 2 

At leading order of perturbation theory in strong coupling constant a s there are 31 
Feynman diagrams contributing to the amplitude of pair J /if) production due to gluon 
fusion. The typical diagrams from this set are presented in Fig. [TJ For the completeness 
we show in Figs. [2H3] also the Feynman diagrams from two other subsets containing 5 and 
8 Feynman diagrams which do not contribute to the production amplitude. Any Feynman 
amplitude shown in Fig. [2]has zero contribution because its color factor is equal to zero. The 
sum of four diagrams from the subset in Fig. [3] is equal zero with the account of relativistic 
corrections studied in this work. So, further we consider only the 31 Feynman amplitudes 
from Fig. [TJ 

Let us consider, for example, the transformation of the first amplitude in Fig. [TJ which 
takes the form: 



r^ h ( Pl , P2 ;q 1 ,q 2 ) = S^ayoe^h^ih 



(f + f+P + ?) (f + f-p-?)' 

u(pih a v{qi)} [u{q 2 )~i v v(p 2 )} 



(7) 
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where £x(ki) and £2(^2) are the polarization vectors of initial gluons. The amplitude ([3]) 
should be convoluted with two wave functions of J/ip mesons taking in the reference frame 
moving with four momenta P and Q. The transformation law of the bound state wave 
function from the rest frame to the moving one with four momentum P was derived in the 
Bethe-Salpeter approach in Ref. 22| and in the quasipotential method in Ref. [23| . We use 
the last one and write the necessary transformation as follows: 



_ d ik p° {R w p)D y* ^«)^(p), 

where R w is the Wigner rotation, Lp is the Lorentz boost from the meson rest frame to a 
moving one, and the rotation matrix D 1;/2 (i?) is defined by 



(8) 



*^(p) 



nl/2/pW 
V l,2 \ K L P 



S-\ Pl , 2 )S(P)S(p) 



where the explicit form for the Lorentz transformation matrix of the four-spinor is 



S(p) 



e(p) + m 
2m 



V«PJ 



e(p) + m 



M 1 "*(k 1 ,k 2 ,P,Q) 



2M5 ab TT 2 al 



dp tt w ( P ) 

(27r) 3 [ e (p) 6(p)+m -i 



J/4, 



(27r) 3 r £ (g) 

I- m 2m J 



(9) 



(10) 



Omitting a number of transformations which can be performed with (!7|) in (131 as in I8L 118 
we write the contribution to the production amplitude in the form: 



Tr< 



v\ — 1 , . p 2 p 

2m(e(p) + m) 2m 



£*p(vi + 1) 



+ Vl 



p^ 



+ 



p 



2m(e(p) + m) 2m 



7°x 



#2 - 1 



+ v 2 



2m(e(q) + m] 



+ 



2m 



$(«2 + 1) 



W 2 + 1 



+ ^2 



xe^(A;i)e§(A;2 



2m(e(q) + m) 
9xv9fi(T 



Q 

2m 
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, 2 • 



where v\ 
and £q • Q 



(f + f + P + g) (f + f-p-g)' 

P/M, t>2 = Q/M; Epq are polarization vectors of final J/"0 mesons with Ep-P 



0; e(p) 



'rm 



+ p 2 , M = 2m is the mass of J/ip meson. The hat 



symbol means contraction of the four-vector with the Dirac gamma-matrices. The spin 
projectors v(0)u(0) = e*(l + 7o)/(2\/2) corresponding to J/tp mesons are introduced as 
well as projectors 5y/v3 onto a color singlet states. We explicitly extracted in ( fill) the 
normalization factors a/2M of quasipotential bound state wave functions. 

The same transformations can be carried out with all 31 Feynman amplitudes in Fig. fT] 
In view of large volume of the calculations we have used the package FeynArts [24j for the 



system Mathematica and Form [25|. To make the entry of final amplitude more compact, 
we introduce a number of vertex functions Tj. Then we can present the total amplitude ([3]) 
in the form: 



M[gg^ J/VJ/V](fci,fc2,P,Q) 



9' 



-MttV 



dp 



dq 



(2tt) 3 J (2tt 



Tr9JT, 



(12) 
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FIG. 2: The subset of LO diagrams which give zero contribution to g + g 
their color factor is equal zero. 



J/ip + J/ip because 






^QQQQQQJ 




FIG. 3: The sum of LO diagrams from this subset gives zero contribution to g + g — > J/ip + J/ip 
with the account of relativistic corrections of order 0(v 2 ). 



[K2 — Pi) — m 



m — ki + pi 
(ki — pi) 2 — m 2 



{K 2 — qi) —m 

r- m + ki — Qi 

(ki — (ft) 2 — m 2 

where inverse denominators of gluon propagators are denned as P>\\ = (k 2 — Pi, 2 — <ft,2) 2 and 
^3,1 = (Pi,2+9i,2) 2 - Pertubative amplitude T(pi,p 2 ; (ft, Q2) in Eq. (02) describes production of 
two c-quark and c-antiquark pairs with the momenta p\^ 2 and q 2 ,i respectively. The formation 
of observable bound states from quark-ant iquark pairs is determined in the quark model 
by the quasipotential wave functions ^ j/^(p, P) and ^j/^(q,Q). These wave functions are 
calculated initially in the meson rest frame and then transformed to the reference frames 
moving with the four- momenta P and Q. As a result we obtain in ffT2]) the following 
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expressions for relativistic bound state wave functions which determine the transition of 
heavy quarks to the bound state: 



*0 J/ *(P) 

e(p) e(p)+m 



L m 2m J 

e p (P,S z )(1 + Vi ) 



Vl + 1 



p- 



2m(e(p) + m 
P 2 

Vl 



V 

2m 



x 



+ 



P 



2m(e(j>) + m) 2m 



(13) 



e(g) e(g)+m 



v 2 - 1 



L m 2m J 

£,(Q,S,)(l + t) 2 ) 



2 

u 2 + 1 



+ #2 



+ 



9 



2m(e(q) + m) 2m 

2 

u 

+ t>2 



X 



9 



'2m(e(g) + m) 2m 



(14) 



Leading order vertex functions in (fT2j) are calculated in the Feynman gauge and can be 
presented as follows: 

r /3 - rn - h + g 2 « _ « m + h - p 2 A 



£2 



(h - g 2 ) 2 - m 2 

m - fc 2 + g 2 
(k 2 - q 2 ) 2 - m 



(k\ — p 2 ) 2 — rn 2 ' 

m + k 2 - P2 
(k 2 - P2) 2 - m 



; £ 2, 



S2 



■ a m + k 2 -p 2 A m - pi - p 2 - gi « 
7 P 77 — —£2 - 8e 27 — — —Y 



(Pi +P2 + gi) 2 - m 2 



m — ki + q 2 
(h - q 2 ) 2 - m 2 L ' (k 2 -p 2 ) 2 -m 

m-k 2 + q 2 r ^ m + k 1 -p 2 m-p l -p 2 -q 1 f 

(k 2 -q 2 ) 2 -m 2 l (ki-p 2 ) 2 -m 2 (pi + p 2 + qi) 2 - m 2 



+ 



(15) 



+ pi + gi + g 2 



m + fci - p 2 „ A m + k 2 -p 2 A 

£27; ^5 n^l + £177 ^5 n £ 2 



18 7a 

V 2 



(Pi + <?i + 92) 2 - "i 2 L (fci - p 2 ) 2 - m 2 (A; 2 - p 2 ) 2 - m 2 

m + fci - p 2 



+ 



m " + " 2 e?7„«S"(pi + g x ) - Pi "'"JV^ i^CPi + *)+ 



(/ci — q 2 ) 2 — m 2 

"' - h • 92 _ a . ;,, 



(k 2 - q 2 ) 2 - m 2 



(fci — p 2 ) 2 — m 2 
m + k 2 -p 2 £ Pn 



(k 2 - p 2 ) 2 - nv 



£2<£i (Pi + 9i ) 



i- 4 — £\ 



m-ki+px r g m + k 2 - qx m - p 2 - gi - g 2 (8 ' 

(&4 — pi) 2 — m 2 L (k 2 — gi) 2 — m 2 (p 2 + gi + g 2 ) 2 — m 2 



+ 



£2 



m - k 2 + p! 
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m + k\ — gi 



•7 



(k 2 — pi) 2 — m? L (fci — gi) 2 — m 2 
g m + pi + p 2 + g 2 



£1 — S£i 



m-p 2 -qi-q 2 ^ 



(P2 + 9i + g2) 2 - "T, 1 



T 



18 7o 



(Pi + P2 + g2) 2 - rn 2 L (A;i - gi) 2 - m 2 
m — ki + pi 



m + k\ — gi . . m + k 2 — qi 

£ 2T1 o £ l + £ 171 ^ o £ 2 



+ 



(fci — pi) 2 — iri 



£?7m^ (P2 + 92) -P2 



(fci — gi) 2 — m 



{k 2 - gi) 2 - m 2 ' 

eMEf (p 2 + g 2 )+ 



8 



m 



fa+Pi 



(fa - Pi) 2 



m' 



m + k 2 - qi 



'(fa 



£ 2 £f a (p 2 + q 2 ) 



q\) L — m z 



l c 2 



s{s v 2 + V 2 r P {pi + qi,p 2 + q 2 ) + V^ v {p 2 + q 2 , Pl + 



m + fa—p 2 
(fa - p 2 ) 2 - m 

m + k 2 -p 2 
X 6 1 (k 2 - p 2 ) 2 - m? 
where we introduce the following tensors: 

1 



a « m-ki + q 2 * 



e 2 -Xe 2 



(fa - q 2 ) 2 - m' 
m-k 2 + q 2 



(fa - q 2 ) 



T P {x, y) = ^xe 1 )(ye 2 )g ap + (fa + x)(k 2 + y)e?4 + e^k" - x a )(2k p 2 - y?)+ 



2x(e 1 e« - e 2 e?)(2A£ - /) -2y{e 1 e{ - e 2 e?)(2A£ 



(16) 



x 



xe x {x a + Ay a )e p 2 - ye 2 (Ax^ + y p )e^. 



Our expressions for the amplitude (|12p and vertex functions (|15p contain relative momenta 
p and q in exact form. In order to take into account relativistic corrections of the second 
order in p and q we expand all inverse denominators of the quark and gluon propagators. 
Such expansions look as follows: 



(Pi + Qi) 2 
1 

m 2 



4 16 

--^ [(p + q) 2 + pQ + qP] +■■ 



4 



(fa - q 2 y 



t - M 2 (t - M 2 f 



[q 2 + 2qk 2 ] + 



(17) 



where the Mandelstam variables for the gluonic subprocess s and t are defined in (jSJ) . There 
are 16 quark and gluon propagators in the amplitude (fT2l) which have to be expanded in the 
same way as in (|T7|) . All denominators of these propagators in nonrelativistic limit take one 
of the following forms: (t — M 2 )/2, (M 2 — s — t)/2, ±s/4 or s/2. Then, the inequalities 



AM 2 < s, 



t + - — M 2 
2 



<£,/i_lM! 

~ 2 V s 



(18) 



mean that in the case of the most unfavorable values of the variables x\ )2 and t we can 
roughly estimate expansion parameters in (TITl) as 2p 2 /M 2 and 2q 2 /M 2 . Preserving in the 
expanded amplitude terms up to the second order in the relative momenta p and q, we can 
perform angular integration using the following relations for iS-wave charmonium: 



*o(p) dp 

r e(p) e(p)+m l ^ 27r) : 
L m 2m J v 7 

*o(p) rf P 



1 



V2tiJ \ 

n L 



p 2 Rs(p) 

fe(p) e(p)+m 



m 2m 



■dp, 



reMe^+ml ( 27r )3 3^ 
L m 2m J 



(fiV - Ul^Vli/) 



p 4 Rs(p) 

e(p) e(p)+m 



m 2m 



dp, 



(19) 



9 



:p-r 
X 



X 



where Rs{p) is the radial charmonium wave function. 

To illustrate the described transformations we present here the result of the calculation 
of the first amplitude in Fig. [1] 

2e P -e* Q (e x -Pe 2 -Q + e x -Qe 2 -P) + 2ep-Q(e x -P e 2 -e* Q + e x -e* Q e 2 -P) - e x -e* P x 

{se 2 -e* Q -2e 2 -Qe* Q -P) - e 2 -s* P (se x -e* Q - 2e x -Qe* Q -P)] (3(1 -c p -c q -c 2 p - c 2 ,) + c p c q x 

(67 + 3c p + 3c q ) + 3c^cy - 64m 2 S [£p-Q(£ r P£ 2 -£Q + £i-£q £ 2 -P) + £q-P(£i-Q£ 2 -£ 

£l-£p£ 2 -Q)] ^3(c p + C q ) + C p C 9 (194 - 3c p - 3c q )j + 16m 2 S 2 [£i-£ p £ 2 -£q + £r£g£ 2 -£ P ] 

(9(c p + c q ) + c p c q (380 - 9c p - 9c ? )) + 192m 2 s£p-4[£ r P£ 2 -Q + e x -Qe r P] (c p + c q + 
c p c q (62 — c p — c 9 )j + 16m 2 s£r£ 2 ^p'^q 32m 2 (3(c p + c q ) + c p c q (329 — 3c p — 3c q )) — 3s 
(3 - 2c p - 2c q - 3c 2 p - 3c 2 q ) + s c p c q (6 13 + 6c p + 6c q ) + 9s c 2 ^] +4£^Q£^P[% + c q ) + 

c p c q (202 - 3c p - 3c q )fj + 512m 2 c p c q ^2 e* P ■ Q e* Q ■ P(e x ■ Q e 2 ■ P + e x ■ P e 2 ■ Q) - 2m 2 e* P ■ e* Q x 

1064m 2 £ r £ 2 + £i-P(125£ 2 -Q - 8£ 2 -P) + £ r Q(125£ 2 -P -8e 2 -Q) 

266 e x ■ e 2 e* P ■ Q e* Q ■ P + e x ■ e* P ( 1 31 s e 2 ■ e* Q + 2 e* ■ P(4 £ 2 • P - 1 29 £ 2 • Q) ) + e x ■ e* Q x 
(I31s£ 2 -£ P + 2£ P -Q(4£ 2 -Q - 129£ 2 -P)) + 2e x -P(Ae 2 -e P e* Q -P - 129 e 2 -e* Q e* P -Q)+ 

2e 1 -Q(4e 2 -e* Q e P -Q - 12$ e r e P e* Q - P)yj^dp dq, 

(20) 

where we introduce the relativistic parameter c p = ^^M- Extracting relativistic factors 

(e+m)/2e in the integrals over both relative momenta p and q we observe that the amplitude 
.M" 6 is a power-like expansion in relativistic parameters c p and c q . Due to the presence of 
four different polarization vectors, which correspond to incoming gluons and outcoming J/ip 
particles, the result (120]) appears to be sufficiently lengthy. We have also obtained analogous 
expressions for remained 30 diagrams, but due to the bulkiness of the total amplitude they 
are not presented here. 

To calculate the cross section we have to sum the squared modulus of the amplitude 
upon all polarizations using the following relations for final J/ip mesons and initial gluons 
correspondingly: 

$>p£r = «-<r, E £ q4^«-^ E<^ = 1 ; 1 k 1 2 -<r- (2i) 

A A A 12 

We find it useful to present the differential cross section for double charmonium produc- 
tion in the proton-proton interaction in the following form: 

^[gg J/ifiJ/iflfrs) = S^l |P(0)| 4 ^a;,P«(t, S ), (22) 

i=0 



m 2 
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where the function describes the LO contribution. It coincides with the nonrelativistic 



analytical expression for the cross section obtained for the studied process in [U|, |13|, Il5|, |26 
The functions F^ 1 ' (i = 1,2,3) describe relativistic corrections. Explicit expressions for all 
functions entering the cross section (|22|) are written in Appendix |A] A number of specific 
parameters uj-i appeared in (122"]) are defined as 



u 



1. 



h 



h 



u 3 



(23) 



They comprise the nonperturbative parameters in the relativistic quark model which deter- 
mine the transition of quarks and antiquarks into the bound states. The parameter -R(O), 
which represents the relativistic generalization of radial wave function at the origin, is defined 
by the formula: 



R(0) 



m + e[p) 

Hp) 



R{p)p 2 dp. 



(24) 



The parameters are determined by integrals containing the bound state wave function in 
the following form: 



m + e(p) 

Hp) 



R(p)p 2 dp, Ji : 



m + e(p) 

Hp) 



m 



e(p) 



m 



e(p) 



1,2 



R(p)p 2 dp. (25) 



Our basic relations for the cross section f)22p evidently show that there exists another 
source of relativistic corrections connected with the charmonium wave functions. For their 
calculation with the desired accuracy we suppose that the dynamics of a cc-pair is deter- 
mined by the QCD generalization of the standard Breit Hamiltonian j27[, which in the cm. 
reference frame can be written as 
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where L = [r x p], S = Si + S 2 , rif is a number of flavors, Ca = 3 and Cp = 4/3 are 
the color factors of the SU(3) color group, je is the Euler constant. The parameter fy of 
vector-exchange confining potential was set to be fy = 0.7. The mass of heavy c-quark in 
our model is equal to m — 1.55 GeV. For the dependence of the QCD coupling constant 
a s (n) on the renormalization point /x in the pure Coulomb term in ( |26i) we use the three- loop 
result [IT 
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whereas in all other terms of the Hamiltonian ([27)1 we use the one-loop approximation for 
the coupling constant a s . The typical momentum transfer scale in a quarkonium is of 
order of the quark mass, so we choose the renormalization scale \i = m — 1.55 GeV and 
A = 0.168 GeV, which gives a s = 0.314 for the charmonium states. The parameters of the 
linear potential A = 0.18 GeV 2 and B = —0.16 GeV have the usual values of quark models. 
Starting with the Hamiltonian ( 1261) we construct the effective potential model based on the 
Schrodinger equation and find its numerical solutions for J/ip meson. Additional details of 
this model are contained in Appendix C of [18j. Note that the obtained charmonium wave 
function is strongly decreasing in the region of relativistic momenta p > m. Our numerical 
evaluation gives the following values of 5-wave charmonium masses: Mj,^ = 3.072 GeV and 
M* h = 2.988 GeV, which lie within the reasonable accuracy in the comparison with their 



L J/i> 
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experimentally measured results Q MfZ = 3.097 GeV and M^ p = 2.980 GeV 



III. NUMERICAL RESULTS AND DISCUSSION 

In this work we investigate the role of relativistic effects in the production of a pair of 
charmonium mesons in proton-proton interaction in the relativistic quark model. We have 
studied only the order a A s parton process of gluon-gluon fusion in the color singlet model. 
At the calculation of the production amplitude (fl2|) we keep relativistic corrections of two 
types. The first type is determined by several functions depending on the relative quark 
momenta p and q arising from the gluon propagators, the quark propagators, and relativistic 
meson wave functions. The second type of corrections originates from the perturbative and 
nonperturbative treatment of the quark-antiquark interaction operator ( 126]) which leads to 
the essential modification of nonrelativistic wave functions. 

For the calculation of relativistic corrections in the bound state wave functions ^ojjl) we 
take the Breit potential (I27p and construct the effective potential model as in 18j, |30| by 
means of the rationalization of the kinetic energy operator. Using the program of numerical 
solution of the Schrodinger equation |3l|] we obtain the following values of all relativistic 

parameters entering the cross section (|2"2"|) : R(0) = 0.57 GeV 2 , U\ = —0.051 and u 2 = 0.0047. 

As it is evident from Eq. (|2"5]) . our definition of integral parameters ii,2 describing rela- 
tivistic contributions from the production amplitude contains the cutoff at relativistic mo- 
mentum of order m. In spite of the convergence of integrals Ii^, our relativistic model can 
not provide a reliable calculation of the wave functions in the region of relativistic momenta 
p > m. So, we introduce a cutoff in (T25|) in order to avoid possible errors caused by the 
mentioned uncertainty. It is obviously that in the quark model we can calculate a number 
of nonperturbative parameters ( 12"3"|) only with certain accuracy. The way of further improve- 
ments in the calculation is related in the first place with more accurate construction of the 
bound state wave function at relativistic momenta. In the approach of NRQCD we encounter 
analogous difficulties connected with the determination of numerous nonperturbative matrix 
elements 0]. 

Let us note also that the cross section (]22|) contains the fourth power of the modified 
wave function at the origin i?(0) and the strong coupling constant a s . Thus, small changes 
of the bound state wave function can lead to substantial changes in final results for the cross 
section. The value -R(O) is calculated with sufficiently high accuracy with the parameters and 
potential ( 12"B"|) of the relativistic quark model. The parameter |i?(0)| 4 undergoes essential de- 
crease in comparison with nonrelativistic value. But other relativistic corrections connected 
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FIG. 4: The differential cross sections for pp — > 2J/tp+X at = 7 TeV as functions of rapidity yp. 
Solid and dashed curves represent total and nonrelativistic results respectively. 



with the functions (i = 1,2,3) in (122]) have the opposite effect on the cross section 
value (|22p . Analytical expression of nonrelativistic contribution to the cross section which is 



determined by F^ coincides with previous calculations in 111. ll3Ml5l. l26j. In the evaluation 
of a s we set the renormalization scale to be the transverse mass /z = rriT = \J ^m? + P^, 
which is generally accepted choice. For the running coupling constant a s (fi) we use the LO 
result with the initial value a s (fi = M z ) = 0.118. 

The basic expression (T5]) for the calculation of the differential cross section contains the 
gluon distribution functions in the proton because the leading contribution comes from 
a gluon fusion process. When the energy of colliding beams increases, the initial parton 
momentum fraction Xj needed to produce heavy quarkonium decreases. It reaches the region 
in x where the number of the gluons becomes much larger than the number of quarks. The 
gluon distribution function determines the probability to find a gluon in the proton with 
some momentum fraction. There exists a number of the parameterizations for partonic 



distribution functions [35|. We use the gluon PDF from the set CTEQ5L as in Ref. [11 



The total numerical value of the cross section obtained from (|2"2"1) is equal to 

a total = 9.6 nb. (29) 

Due to law-x behavior of gluon distribution functions, the main contribution to the integral 
cross section (|2"9"1) results from the region x\ 2 ~ 10~ 3 . To be more precise, this region 
is determined by the condition: 7.8 ■ 10~ 7 < X1X2 < 7.8 ■ 10 -6 . In the frequently used 

nonrelativistic limit when i/;(Q) = 0.21 GeV^ and R(0) = 0.74 GeV 5 the total value of 
the cross section amounts a nonre i = 18.3 nb and agrees with the calculation in (TTj | . In the 
nonrelativistic limit of our quark model based on Eqs. f f26|) and f l27|) we obtain slightly 

_ 3 

greater values -R(O) = 0.79 GeV 2 and a nonre i = 23.1 nb. To obtain ( 1291) the factorization 
scale in the parton distribution functions f g / p (x,fi) is taken equal to the transverse mass 
too: fi = uit- 

To compare the results of our calculation with the measured value of the cross section 

in 



10J it is necessary to write the differential cross section in terms of the rapidity yp 



\ In jrzrpj ■ The rapidities of outcoming charmonia with momenta P and Q can be obtained 
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Pr, GeV F T , GeV 

FIG. 5: The differential cross sections for pp -> 2 J/ip + X at ^ = 7 TeV (left) and VS = 14 TeV 
(right) as functions of transverse momentum of the J/i/? pair integrated over the rapidity. Solid 
and dashed curves represent total and nonrelativistic results respectively. 



in the form: 

1 1 Xl JL 1 1 

^ = 2 ln ^ ± 2 ln 



M 2 -t 



(30) 



The differential cross section da/dyp for the reaction pp — > 2J/tp + X is shown in Fig. HI 
It is clear from this plot that relativistic effects strongly influence on the rapidity distribu- 
tion of the final charmonium. In the LHCb experiment [Io| the rapidity lies in the range 
2 < Up,q < 4.5, so we should integrate the differential cross section fl2]) over rapidities from 
such interval in order to obtain the value corresponding to the experiment at the LHCb 
detector. Then we obtain: 

a th {2 < y P;Q < 4.5) = 1.6 nb. (31) 

The result (13 ip is significantly smaller than the experimental value of the cross section ([T]). 
Different sources of relativistic corrections in ([2]) are differently directed. But we observe 
that combined action of all relativistic effects lead to essential decreasing of the production 
cross section. In this work we carry out the investigation only of one important source 
of corrections to the nonrelativistic cross section. Decreasing behavior of the cross section 
a(pp — > 2J/ip+X) due to the account of relativistic contributions is noticeable clearly in spite 
of existing theoretical errors occurred in our calculation. In our analysis of the production 
amplitudes we correctly take into account relativistic contributions of order 0(v 2 ). Therefore 
the first basic theoretical uncertainty of our calculation is related with the omitted terms 
of order 0(t> 4 ). Since the calculation of the charmonium mass is sufficiently accurate in 
our model (the error is less then 1%), we suppose that the uncertainty in the cross section 
calculation due to omitted relativistic corrections of order 0(v 4 ) in the quark interaction 
operator (the Breit Hamiltonian) is also very small. Taking into account that the average 
value of the heavy quark velocity squared in the charmonium is (v 2 ) = 0.3, we expect that 
relativistic corrections of order 0(v 4 ) to the cross section (!3T|) coming from the production 
amplitude should not exceed 30% of the obtained relativistic result. As we mentioned above 
in the quasipotential approach we can not find precisely the bound state wave functions in 
the region of relativistic momenta p > m. Using indirect arguments related with the mass 
spectrum calculation we estimate in 10% the uncertainty in the wave function determination. 
Larger value of the error will lead to the essential discrepancy between the experiment and 
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theory in the calculation of the charmonium mass spectrum. Then the corresponding error 
in the cross section (13T|) is not exceeding 20%. We do not consider a part of theoretical error 
related with radiative corrections of order a s because these corrections are omitted in our 
analysis. So, our total theoretical error is not exceeding 36%. To obtain this estimate we 
add the above mentioned uncertainties in quadrature. 

We show in Fig. [5] the distribution over transverse momentum of the J/ip mesons in- 
tegrated over all rapidities at v^S* = 7 TeV (left) and yS = 14 TeV (right). Previous 
investigations [HI, EH, 14] of the pair charmonium production in pp interaction showed that 
the color singlet channel prevails in the differential cross section do jdP^{j>p — > 2 J ftp + X) 
at small Pt, but the color octet channel dominates at large Pt- It can be seen in Fig. [5] 
at \/~S = 7 TeV that the account of relativistic corrections leads to the ratio of relativistic 
and nonrelativistic cross sections cr re i/a nr ~ 0.4 near the peak. This trend remains un- 
changed in the region of high transverse momenta. So, the color-octet contribution retains 
the dominance at large Pt- We investigate also the relative value of relativistic corrections 
in the production rate with the growth of the energy v^S*- Our calculation show that at 
VS = 14 TeV (see the right plot in Fig. |5]) the ratio of relativistic and nonrelativistic cross 
sections is retained without essential modifications. The cross section increases with the 
growth of the energy and reaches the value a th (2 < y Py q < 4.5) = 2.98 nb at \/S = 14 TeV. 
It is appropriate to mention here one result regarding to the study of relativistic effects in 
single J/ip production at hadron collisions j32|. It was shown in that paper that relativistic 
corrections to the color-singlet J/ip hadropro duct ion of order 0(v 2 ) are at a level of about 
1% for sufficiently large Pt'- 5 < Pt < 50 GeV. Our calculation demonstrates that in the 
region of transverse momenta 5 < Pt < 50 GeV the value of relativistic corrections in 
the cross section of the pair charmonium production reaches 60%. Relativistic corrections 
which we study in this work include not only the terms of order 0(v 2 ) in the production 
amplitude but also the same order effects in the long-distance matrix elements. In spite 
of the difference between ([T]) and (j3"Tl) . we consider that at present it is difficult to state 
that there is the discrepancy between the theory and experiment in double charmonium 
production. Indeed, it is known that NLO in a s contributions have large value in inclu- 
sive single- quarkonium production at hadron colliders 0, IH, 34]- The example is found in 
the inclusive J/ip production where the NLO corrections to the color-singlet contribution 
increase the total cross section by a factor of about 2 and the production rate of J/ip is 
much increased for larger transverse momentum Pt- Therefore, one can expect that the 
NLO corrections to the double charmonium production in proton-proton interaction can 
smooth the appeared difference between (pQ) and (13"Tj) . Moreover, as we mentioned above 
there exists new mechanism through the double parton scattering which gives the contribu- 



tion comparable with the standard nonrelativistic result: (Jdps(pp — > 2J/vp + X) = 2 nb [16 
Accounting for this result and our value of the cross section f )3T|) we obtain the summary 
value a(pp — > 2J/tp + X) = 3.6 nb. Then, taking into account the experimental error, the 
difference with the LHCb experiment does not look so significant. 
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Appendix A: The coefficients entering the differential cross section (1221) 
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2073684017664s 3 * 13 + 987205294080s 2 * 14 + 288123568128s * 15 + 38578913280* 16 ) . 
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